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The problem of error analysis is addressed in stages beginning with the case of uncorrelated 
parameters and proceeding to the Bayesian problem that takes into account all possible correla- 
tions when a great deal of prior information about the accessible parametr space is available. The 
formulas for the standard deviations and deviations with arbitrary confidence levels are derived. 
Underestimation of the errors of XAFS-function extraction is shown to be a source of unjustified 
optimistic errors of fitting parameters. The applications of statistical \ 2 ~ an< A f-tests to the fitting 
problems are also discussed. 
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£ , I. INTRODUCTION 

In the Open Letter to the XAFS Community [1] Young and Dent, the leaders of the UK XAFS User Group, 



expressed their concern over the persistence of lingering common opinion that XAFS is a "sporting technique" and it 
is possible to obtain the "answer you want" . Some way out they see in a special attention to the publishing XAFS data 
(first of all, to XAFS spectra) and have formulated several recommendations for editors and referees. Undoubtedly, 
, in the matter of extraction of the real, not invented, information from XAFS experiments the quality of spectra is 
of great importance. We see here another problem as well. Not having some necessary elements of XAFS analysis 
£>^, (some values and the procedures for their determination), one has a quite natural desire to turn those values to 
advantage. Principally we mean the inability of the standard methods to find the errors of the atomic-like background 
C^j /io- Traditionally, the noise is assigned to these errors. However, as was shown in Ref. [2], the noise is essentially lower 
than the errors of the /io construction. Below, we will show that the underestimation of the errors of XAFS-function 
,_H ' extraction is a source of the unreasonable optimistic errors of fitting parameters. 

. Practically all known programs for XAFS modeling [3] in some way calculate confidence limits of fitting parameters. 
However, since there is no standardized technique for that and since most published XAFS works do not contain any 
mention of methods for estimation of the errors of fitting parameters, the accuracy of the XAFS results remains to 
be field for trickery. 

In the present article we derive the expressions for the errors of fitting parameters under different assumptions on 
the degree of their correlation. Besides, the prior information about parameters is possible to take into account in the 
framework of Bayesian approach. Moreover one can find the most probable weight of the prior information relative 
""^ . to the experimental information. 

We also discuss the grounds and usage of the statistical tests. The special attention was focused on that where and 
* ^ , how one can embellish the results and artificially facilitate the statistical tests to be passed. 
J^sj' All methods and tests described in the paper are realized in the program VIPER [6]. 



II. ERRORS IN DETERMINATION OF FITTING PARAMETERS 



Let for the experimental curve d defined on the mesh x\,. . . ,xm there exists a model m that depends on N 
■ parameters p. In XAFS fitting problems as d may serve both (not weighted by k w ) and x( r )- The problem is 
to find the parameter vector p that gives the best coincidence of the experimental and model curves. Introduce the 
figure of merit, the x 2 -statistics (do not confuse with the symbol of XAFS function): 

i=i £ i 

where £$ is the error of di. The variate x 2 obeys the x 2 -distribution law with M— N degrees of freedom. Of course, for 
the given spectrum d and the given model m the value of \ 2 is fully determined; we call it "variate" bearing in mind 
its possible dispersion under different possible realizations of the noise and the experimental errors of di extraction. 
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Often a preliminary processing (before fitting) is needed: smoothing, filtration etc. Naturally, during the pre- 
processing some part of the experimental information is lost, and on the variates = (di — m^jsi additional 
dependencies are imposed (before, they were bound solely by the model m). It is necessary to determine the number 
of independent experimental points iV; n d. For the commonly used in XAFS spectroscopy Fourier filtering technique 
the number of independent points is given by [4]: 

N ind = 2AkAr/ir + 2, (2) 

where Ak = fc max — fc m ; n and Ar = r max — r m ; n are the ranges in k- r-spaces used for the analysis, and r m ; n > 0. If 
r m in = then 

N ind = 2AfcAr/7r + f . (3) 

Instead of keeping in the sum (1) only A ind items which are equidistantly spaced on the grid X\, . . . ,x M , it is more 
convenient to introduce the scale factor N- ln d/M: 

2 ^ind (<k - mi) 2 

x= t^^t- (4) 

i=l % 

Now the variate \ 2 follows the ^-distribution with N[ n< i — N degrees of freedom. It can be easily verified that with 
the use of all available data (r m ; n = and r max = n/2dk) the definition (4) turns into (I). 

Let us now derive the expression for the posterior distribution for an arbitrary fitting parameter pj : 

P( Pj \d) = J •••<//>,,,•••/': P<F. (5) 

where P(p|d) is the joint probability density function for all values p, and the integration is done over all Pijtj. 
According to Bayes theorem, 

P(d|p)P(p) 

P (Pl d ) = pTj) ' ( 6 ) 

P(p) being the joint prior probability for all pi, P(d) is a normalization constant. Assuming that -/V; nd values in d 
are independent and normally distributed with zero expected values and the standard deviations £j, the probability 
P(d|t), so-called likelihood function, is given by 

P(d|p) cx exp (- X 2 /2) , (7) 

where x 2 was defined above by (4). Its expansion in terms of p near the minimum (V p % 2 = 0) which is reached at 
P = P yilds: 

i i N pa 2 

P(d|p) cx exp(- i( p - p) T ■ H ■ (p — p)) ^ (- i ^ ^Ap k Ap) , (8) 



where Apk = Pk — Pk, and the Hessian H components (the second derivatives) are calculated in the fitting program 

1 kl 



at the minimum of x 2 . The sufficient conditions for the minimum are Hkk > and H/j/jHh — H 2 ; > 0, for any k,l. 



Hence, the surfaces of constant level of P(d|p) are ellipsoids. 



A. Simplest cases 

If one ignores the prior then the posterior probability density function P(p|d) coincides with the likelihood P(d|p). 
Let us consider here two widely used approaches. 

(a) Parameters are perfectly uncorrelated. In this case the Hessian is diagonal and 

P(p,|d)cxcxp(-iH JJ Ap 2 ). (9) 

The standard deviation of pj is just 
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6^ Pj = (2/H w ) 1/2 . (10) 
(b) Parameter pj essentially correlates solely with pi. In this case 

P(pj\d) oc J dpiPfapjld) oc J d Pi exp(-^ii jj (A Pj ) 2 - ^H^ApjAp, - |H^(Ap 4 ) 2 ) 

ex exp(-l [h# - SyS i4 ] (Ap,) 2 ) , (11) 
from where one finds pj = pj and the mean-square deviation 



In practice, to find the strongly correlated pairs of parameters, one finds the pair-correlation coefficients: 

(13) 



(ApiApj) - (A Pi )(A Pj ) 



6(A Pi )6(A Pj ) 

taking on the values from -1 to 1. Two parameters are uncorrelated if their correlation coefficient is close to zero. It 

is easy to calculate the average values over the distribution (11): (Ap 2 ) = 2Hjj/Det, {Ap 2 ) = 2Hj,/Det, (ApiApj) = 

— 2Hjj/Det, where Det = YijjUa — H 2 ^. Notice, by the way, that these arc the elements of the inverse matrix of H/2. 
Now the pair-correlation coefficients are given by: 

- (») 



Via the correlation coefficient the mean-square deviations, found for the cases (a) and (b), are simply related: 



6^ Pj =6^p J\-r%. (15) 



Consider an example of the error analysis. For L 3 Pb absorption spectrum 1 for BaPbC-3 compound the average error 
of the XAFS extraction from the measured absorption was £j = 0.007. For the filtered over the range 1.0 < r < 2.1 A 
(the signal from the octahedral oxygen environment of lead atom) XAFS (see Fig. 1), the model function was calculated 
as follows. For one-dimensional the Hamiltonian of the lead-oxygen atomic pair with potential U = a/2 ■ (r — ro) 2 we 
found the energy levels and corresponding to them wave functions. Then, averaging over the Gibbs distribution, the 
pair radial distribution function (PRDF) normalized to the coordination number TV was found as: 

g(r) = Nj2\*n(r)\ 2 e- E -/ kT /j2e- E " /kT , N= j g(r) dr, (16) 

n n J 

and the XAFS function as: 

^max 

x (k) = ^F(k) J g(r)sm[2kr + cj)(k)}/r 2 dr. (17) 

f min 

The phase shift 4>(k) and the scattering amplitude F(k) were calculated using feff6 program [5]. By variation of 
the parameters ro, a, TV (where TV includes the factor S 2 ), and Eq, the shift of the origin for the wave number k, one 
search for the best accordance between the model and experimental curves. Here for the fitting, the VIPER program 
was used which, in particular, calculates the Hessian of x 2 (defind by (4) with TVi n d = 11.8) at the minimum. The 
correlation coefficients are listed in the Table I. 



1 The spectrum was recorded at 50 K in transmission mode at D-21 line (XAS-13) of DCI (LURE,Orsay, France) at positron 
beam energy 1.85 GeV and the average current ~ 250mA. Energy step — 2eV, counting time — Is. Energy resolution 
of the double-crystal Si [311] monochromator (detuned to reject 50% of the incident signal in order to minimize harmonic 
contamination) with a 0.4 mm slit was about 2-3 eV at 13 keV. 
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TABLE I. Pair-correlation coefficients rtj for the example fitting. 

N a ro Eo 

~N I -0.286 O092 0M1 

a -0.286 1 -0.044 0.048 

r 0.092 -0.044 f 0.727 

E 0.041 0.048 0.727 1 



TABLE II. Mean values and mean-square deviations of the fitting parameters. Sp are the mean-square deviations calculated: 
for perfectly uncorrelated parameters (a), trough the maximum pair correlations (b), from the bayesian technique without prior 
information (maximum likelihood) (c), from the posterior probability that considers the most probable contribution of the prior 
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We now turn our attention to the errors of fitting parameters. In ignoring the correlations, the errors S^p are rather 
small (see Table II). However, we know that the parameters r and E are highly correlated, and their real errors 
must be larger. In the traditional XAFS-analysis two-dimensional contour maps have long been used [7] for estimates 
of the correlation score and the error bars. Notice, that to do this requires, first, the definition and determination of 
the correct statistical function x 2 (but not a proportionate to it), and, second, a criterion to choose the critical value 
of x 2 (depending on the chosen confidence level). 

For the most correlated pair, ro and Eq, find the joint probability density function P(ro-Eo|d) using the Hessian 
elements found at the minimum of the \ 2 : 

P(r E \d) cx cxp(-iH roro (Ar ) 2 - ^H roEo Ar AE - ^H Eo e (A£ ) 2 ) (18) 
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FIG. 2. The joint probability density function P(roEo\d) calculated via the expansion (8) (solid lines) and using the exact 
X 2 function (on the right, dashed lines). Also shown the graphical interpretation of the mean-square deviations <^ a Vo and 
ro given by (10) and (12). The ellipse of the standard deviation is drawn by the thick line. 

which is depicted in Fig. 2 as a surface graph and as a contour map. The ellipses of the equal probability are described 
by: 

H roro (Ar ) 2 + 2R roEo Ar AE + R EoEo (AE Q ) 2 = 4A. (19) 

In Fig. 2 they limit such areas that the probability for the random vector (ro,Eo) to find itself in them is equal to 
£ = 1 — e~ x = 0.2, 0.6, 0.8, 0.9 and 0.95. By the thick line is drawn the ellipse corresponding to the standard deviation: 
A = 1/2 and 1=1 — e -1 / 2 « 0.3935. For this ellipse the point of intersection with the line AE = and the point 
of maximum distance from the line Aro = give the standard mean-square deviations S^tq and S^tq that coincide 
with the expressions (10) and (12). To find the mean-square deviation 5^' for an arbitrary confidence level t, one 
should multiply the standard deviation by — 21n(l — £). 

In Table II the errors in the column S^p were found as the largest errors among all those calculated from the pair 
correlations. For the parameters N and a all pair correlations are weak, so their 5^ and <5^ b ^ are hardly differ. For 
the parameters ro and Eq these mean-square deviations differ remarkable. 

Finally, we put the question, how much is rightful the expansion (8) for the likelihood function? In Fig. 2, on the 
right, the dashed ellipses of equal probability are found for the exact \ 2 that was calculated by the viper program 
as well. Mainly, just-noticeable difference is caused by the realization of the fitting algorithm or to be more precise, 
by the values of the variations of the fitting parameters which determine the accuracy of the minimum itself and the 
accuracy of the derivatives at the minimum. Of course, this difference can be neglected. 

B. General case 

Often, a particular fitting parameter significantly correlates not with a one, but with several other parameters 
(in our example this is not so, but, for instance, the problem of approximation of the atomic-like background by 
interpolation spline drawn through the varied knots [8,2] is that very case). Now, the consideration of the two- 
dimensional probability density functions is not correct no more, one should search for the total joint posterior 
probability P(p|d). 

For that, first of all, one is to find the prior probability P(p). Let we approximately know in advance the size Sk 
of the variation range of the parameter pk- Then the prior probability can be expressed as: 

N A 2 

P(p| a )«^/ 2 exp(-f]r^), (20) 

k=l k 
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where the regularizer a specifies the relative weight of the prior probability; at a — there is no prior information, 
at a — > oo the fitting procedure gives nothing and the posterior distribution coincides with the prior one. In the 
expression (20) a appears as a known value. Later, we apply the rules of probability theory to remove it from the 
problem. 

So, for the sought probability density functions we have: 

f 1 N 

P(pj\d, a)ct ■■■ dpi^j ■ ■ ■ a N/2 expf-- g kl Ap k Ap^ , (21) 

J k,l=i 

where 

a £ , Hfei , . 

9ki = ^2<3fci + — ■ (22) 

Since there is no integral over pj, separate it from the other integration variables: 

N N 



1 f 1 

P(Pj\d, a) cx a N/2 exp(--gjjApfj ■■■ dp^ ■ ■ ■ exp(-- ^ g k iA Pk A P i - A Pj ^ g kj Ap^J , 

u l i i 1 



(23) 



fcj=l k=l 



Here, the symbol j near the summation signs denotes the absence of j-th item. Further, find the eigenvalues Aj 
and corresponding eigenvectors e, of the matrix g k i in which the j-th row and column are deleted, and change the 
variables: 

N N , 

h = y/\iYf Ap k e ik , Ap k =Yf^ (i,k^j). (24) 
k=i <=i VAj 

Using the properties of eigenvectors: 

N N 

J~] J gikZik = Xjeu, e lk e ik = 5 H (25) 

k=l fe=l 

one obtains: 

P( Pj \d, a) cx a N ' 2 cxp(--[ fe - v?]Apf) ■■■ db w ■ ■ ■ exp(-- ^% + Ui A Pj } 2 ) 

cxa N / 2 exp(-i[. fe -u 2 ]A P 2 ), (26) 
where new quantities were introduced: 



N N 



Ui = -^=YJ g kj e ik , u 2 ^"7- (27) 



Thus, we have found the explicit expression for the posterior distribution of an arbitrary fitting parameter. This is 
a Gaussian distribution with the mean pj = pj and the standard deviation 

S^ Pj = ( gjj -u 2 )- 1/2 - (28) 

The formulas (26)-(28) require to find the eigenvalues and eigenvectors for the matrix of rank N — 1 for each 
parameter. Those formulas have merely a methodological value: the explicit expressions for posterior probabilities 
enables one to find the average of arbitrary function oipj. However, the standard deviations could be calculated 
significantly easier, having found the eigenvalues and eigenvectors for the matrix of rank N one time. 

(c) _ /A^Pfald, _ J Ap 2 cxp(-lEM = iff^A Pfc Ap,)rfp 

Pj " lP( Pj \d,a)d Pj ~ /exp(-|EM=iff«Ap fc A Pl )dp 
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Analogously to what was done above, performing the diagonalization of gki, one obtains: 



/dbexp(-jE^i^) 



X, ' 



i=l 



(30) 



where the eigenvalues (\) and eigenvectors (ej) correspond to the full matrix g^x- 

One can give another interpretation of the S^p-Hnding process. It is easy to verify that H/2 and the covariance 
matrix C of the vector p are mutually inverse. Therefore 



(8 {c) Pj) 2 = C jj =2(K- 1 ) jj , 



(31) 



and the variate (p — p) T • C _1 • (p — p) = |(p — p) T ■ H • (p — p) is x 2 -distributed with N degrees of freedom if p is 
the TV-dimensional normally distributed vector (by Eq. (26) this condition is met). The ellipsoid that determines the 
standard deviation is: 



(p - p) T • H • (p - p) 



N. 



(32) 



For an arbitrary confidence level £, on the r.h.s. would be {xn)£, the critical value of the x 2 -distribution with N 
degrees of freedom. The error d^pk is equal to the half the ellipsoid size along the fc-th axis. 

In our example fitting, the errors found in the absence of any prior information (a = 0) from the formula (30) are 
listed in Table II in the column S^p. Due to every one parameter correlates at the most with one other parameter, 
all S^p are practically coincide with 5^p. Generally, this may be not so. 

Finally, let us find the most probable value of a. Its posterior distribution is given by: 



P(a|d) = J d P P(a,p\d) = J dpP(a)P(p\a,d). 
Using a Jeffreys prior P{a) — 1/a [9], one obtains for the posterior distribution: 

f 1 N 

P(a\d) oc / dpa^/^exp^-- 9ki^Pk^Pi) oc (Ai^-Ajv)" 



l/2 a JV/2-l_ 



(33) 



(34) 



k,i=i 



In our example we have set the variation range of the parameter pk to be equal to Sk — ±Pfc (this means that 
Pk G [0, 2pk\) for all parameters except for Eq', since it varies near zero, we have chosen Se = ±10 cV. For the 
mentioned variation ranges, the distribution P(a\d) has its mode at a = 2.64 • 10 3 (see Fig. 3). The bayesian errors 
found for this regularizer are listed in the column <5^ d ^ of Table II. As a result, we have got the mean-square errors 
that for some parameters are significantly lower than even S^p. There is nothing surprising in that: any additional 
information narrows the posterior distribution. If we would choose Sk to be less, S^'pk would be yet lower. For 
instance, XAFS is quite accurate in distance determination, and for many cases one can assume distances to be 
known within ±0.2 A. In our case this leads to 6^r = 3.4 • 10~ 3 A. 




FIG. 3. The posterior distribution for the regularizer 
a found from Eq. (34). 
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C. Important note 

Having obtained the expressions (10), (12) and (30) for the errors of fitting parameters, we are able now to draw 
an important conclusion. If in the definition (4) one substitutes for £j the values that are smaller by a factor of (3 
than the real ones, the x 2 and its Hessian's elements are exaggerated by a factor of /3 2 , and from (10), (12) and (30) 
follows that the errors of fitting parameters are understated by a factor of (3\ 

In the preceding paper [2] it was shown that the errors of the atomic-like absorption construction are essentially 
larger than the experimental noise, and therefore it is the former that should determine the £j values. However, 
these values are traditionally assumed to be equal to the noise, or one uses unjustified approximations for them, also 
understated (like 1/ef = k w [10]). It is here where we see the main source of the groundless optimistic errors. 

III. STATISTICAL TESTS IN FITTING PROBLEMS 

A. x 2 -test 

Introducing the statistical function x 2 , we assumed that it follows the x 2 distribution with v = M — N degrees of 
freedom. However for this would be really so, one should achieve a sufficient fitting quality. This "sufficient quality" 
could be defined as such that the variate (4) obeys the x 2 distribution law, that is this variate does not fall within 
the tail of this distribution. Strictly speaking, the following condition must be met: 

X 2 < (xl)e, (35) 

where the critical value (xl)e f° r the specified significance level I may be calculated exactly (for even v) or approxi- 
mately (for odd v) using the known formulas [11]. 

Notice, that the choice of the true Si here also plays a cardinal role. However, it is important here that one would 
not use the overestimated values which facilitate to meet the requirement (35). As we have shown in [2], one could 
obtain the overestimated £j, having assumed the Poisson destribution law for the detectors counts when the actual 
association between the probability of a single count event and the radiation intensity is unknown. 

Thus, the exaggerated values tell about a quality fitting, but give the large errors of fitting parameters. The 
understated Si lead to the would-be small errors, but make difficult to pass the x 2 -test (i. e. to meet the condition 
(35)). We are aware of many works the authors of which do not describe explicitly the evaluation process for the 
errors of XAFS-function extraction and do not report their explicit values. However, by implication it is seen that Si 
were chosen (not calculated!) as low as possible to scarcely (with t = 0.9 — 0.95) pass the x 2 -test; as a result, very 
impressive errors of the structural parameters were obtained. In such approach no wander that the difference of 0.01 A 
between the diffraction data and the XAFS-result that was found within 0.002 A was attributed to the "suggested 
presence of a small systematic error" [10]. 

B. F-test 

Let there is a possibility to choose between two physical models depending on different numbers of parameters Aq 
and N2 (N2 > Ni). Which one of them is more statistically important? For instance one wish to decide whether a 
single coordination sphere is split into two. 

Let for the two models the functions x? and xl obey the x 2 -distribution law with v\ = Ni n( ± — Aq and v 2 — N- m d — N 2 
degrees of freedom, correspondingly. From the linear regression problem (near the minimum of x 2 > the likelihood 
function is expressed by (8) and is identical in form to that of the linear regression problem) it is known that the 
value 

, (xl - xDliyi - vn) 

1 Xll»2 

obeys the Fisher's F-distribution law with (y\ — V2, v-i) degrees of freedom if exactly r — v\ — v 2 parameters in the 
second model are linearly dependent, that is if exist the rxJV 2 matrix C of rank r and the vector c of the dimension 
r such that Cp = c. In order for the linear restrictions on the second model parameters to be absent, the value / 
should not follow the ^-distribution, that is it should be greater than the critical value {Fv 1 -v 2 ,v 2 )i f° r the specified 
significance level i: 
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(36) 




/ > (F Vl _ V2i „ 2 )e (37) 

or 

xl < Xl ({F^-^)"- 1 ^ 1 + lj . (38) 

Notice, that the expression (38) means the absence of exactly r linear restrictions on the second model parameters. 
Even if (38) is realized, the less number of linear dependencies are possible. If, for instance, the splitting of a single 
coordination sphere into two does not contradict to the F-test (38), some of the parameters of these two spheres may 
be dependent, but not all. This justifies the introduction of a new sphere into the model XAFS function. 

Thus, having specified the significance level £, one can answer the question "what decrease of x 2 must be achieved 
to increase the number of parameters from Ni to A 2 ?" or, inside out, "what is the probability that the model 2 is 
better than the model 1 at specified (Ni,xt) an d (A^xi)?" 

Notice, that since in the definition for / the ratio x\lx\ appears, the actual values of £j become not important for 
the F-test (only if they all are taken equal to a single value). 

Consider an example of the statistical tests in the fitting problem. In Fig. 4 are shown the experimental curve 
with iVi n d = 11.8 and two model curves with Ni = A and N2 — 7. The underlying physical models were described in 
Ref. [12]; here only the number of parameters is of importance. Let us apply the statistical tests. Through the fitting 
procedure for the model 1 we have: v\ = 11 — 4 = 7, Xi = 16.8 > 14.1 = (x?-)o.95, for the model 2: = 11 — 7 = 4, 
Xi = 5.3 < 9.5 = (x!)o.95- That is the first model does not pass the x 2_ t es t- Further, / = 2.89 = (^3,4)0.84, from 
where with the probability of 84% we can assert that the model 2 is better than the model 1. 

In the XAFS analysis the F-test has long been in use [7]. However, the words substantiating the test are often 
wrong. The authors of Rcfs. [10,13], for example, even claimed that the value / (36) must follow the F-distribution, 
although then in Ref. [13] there appears a really correct inequality (38). 



IV. CONCLUSION 



The solution of the main task of the XAFS spectroscopy, determination of the structural parameters, becomes 
worthless if the confidence in this solution is unknown. Here we mean not only the confidence in the obtained fitting 
parameters that is their mean-square deviations, but also the credence to the very methods of the error analysis. It 
is excessive optimistic errors evaluations lead to the suspicious attitude to the XAFS results. 

To improve the situation could the development of the reliable and well-grounded techniques that do not allow 
one to treat the data in an arbitrary way. First of all, this is a technique for determination of the real errors of the 
atomic-like absorption construction. Second, we regard as necessary to standardize the method for the correct taking 
into account of all pair correlation between fitting parameters. And third, (we have not raised this question here) 
programs for scattering phase and amplitude calculations should report on the confidence limits for the calculated 
values, that is report how sensitive the calculated values are to the choice of the parameters of scattering potentials. 
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